%% OSLWORKSHOPS - DEMONSTRATING HOW TO CLEANUP DATA USING OSLVIEW & AFRICA

% This is a practical to intorduce you to OSLVIEW and OSL_AFRICA.
%
% OSLVIEW is a MEG data viewer for continuous data.
% OSL_AFRICA is OSL's de-noising program for MEG data.

% In this practical you will be cleaning up 5 minutes of resting state data
% acquired on an Elekta Neuromag. Maxfilter has been applied to the data.

% You will

% 1.) Use OSLVIEW to flag bad channels and bad periods of data.

% 2.) Use OSL_AFRICA to automatically detect the EOG, ECG and MAINS
% artefacts. You will also manually inspect the 5 most and 5 least kurtotic
% components. 

% 3.) You will compare the orginal and denoised data to assess the impact
% of the denoising.

% You will need OSL version 1.2.beta.2 or later. 

% Henry Luckhoo
% henry.luckhoo@trinity.ox.ac.uk
% 06.11.12




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

error('THIS EXAMPLE HAS NOT BEEN UPDATED TO WORK WITH OSL1.3.1');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%% Setup Paths for your machine - Edit and Run this cell.

tilde='/home/gileslc';
osldir=[tilde '/OSL-Repo/osl/']; % Change this to be the path to your version of OSL.    
rmpath(genpath('/net/horus/data/OHBA/spm8')); % If you have other versions of SPM on your path you should remove them too!
addpath(osldir);
global OSLDIR;
osl_startup(osldir);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set-up files and directories - Edit and Run this cell.

workingdir=['/home/scratch_disk/hluckhoo/OSL_WORKSHOP_AFRICA_practical_data']; %  % Change this to be the path to the practical data directory.
cd(workingdir);

spm_files = {[workingdir '/spm8_rest_data_subj01_sss']};                % This is a list of the SPM files you plan to denoise.
eye_tracker_files = {[workingdir '/rest_data_subj01_eyetracker.txt']};  % These are the corresponding Eyelink Eyetracker files.

%% A note on EyeLink Eyetracker Data.

% If you have acquired eyetracker data, you can use it to identify blink
% related activity. If you have recorded Eyetracker (both monocular and
% binocular will work), an .edf file will be stored on the EyeLink
% computer. You need to copy this data off the EyeLink computer, on to the
% stimmeg computer (via the shared network folder). On stimmeg, use the 
% EyeLink software to produce a report file. I can show you how on request.
% This gives a txt file with the necessary data for blink detection. In
% principle this could work for saccades as well but I haven't attempted
% it.

%% Visual Inspection of the Data - Run the cell and follow instructions.

% 1.) Run this cell to load an SPM object into OSLVIEW. 
subnum = 1;             
D = spm_eeg_load(spm_files{subnum});
oslview(D);

% 2.) In the GUI, change the time axis to look at 2 or 3 minutes of data by
% clicking on the "Increase time window button" in the top left (note that
% the x-axis is labelled in seconds). 

% 3.) Using the variance bars on the left, set any noisy planar
% gradiometers to bad by right clicking on the noisy channel and selecting 
% "Set Channel as Bad". There are two clear outlier channels. Set them to "bad".

% 4.) You can alway right click on a bad channel and set it to "good" by 
% clicking on "Set channel as Good". Try setting a different channel to
% "bad" and then back to "good".

% 5.) There is an obvious artefact in the data. Find it and mark that
% period as "bad". To mark noisy periods as bad epochs you need to right
% click before the bad data and selecting "Mark Event". Do the same to mark
% the end of the bad period. You must always mark the start and end of a
% bad epoch. As soon as the you mark the end of the bad epoch, the epoch
% will be marked as bad, and the corresponding data disappears from the
% average channel at the bottom.   

% 6.) Note that if you want to undo the marking of a bad epoch, you can
% right clicking on the data between the start and end markers and
% selecting "Remove Event". 

% 7.) Click on the "P" button at the top of the GUI. This switches the
% sensor types to the magnetometers "M". Mark any bad magnetometers as bad.
% Also mark any bad epochs that weren't seen in planar gradiometers.

% 8.) Click on the save button at the top. Then close the GUI. When
% prompted to save the figure, click "No".

%% Look at the updated SPM object

% If you have used OSLVIEW correctly there should be some bad channels and
% bad epochs saved in the SPM object. Run this cell to check that you have
% done this successfully.

clear D
D = spm_eeg_load(spm_files{subnum});

% load in the events
eve=events(D);

% to see event 1:
eve(1)

% to see the bad channels
D.badchannels

%% Now let's use AfRICA to clean up the data.

% AfRICA has 3 parts.

% 1.) ICA decomposition of the MEG data.
%     - Use recommended settings.
%     - You should specify that it is using Maxfiltered data.
%     - This stage can be slow but can be run over multiple subjects and
%       saved to disk (e.g. overnight, over weekend).
%     - To do this, set all stage 2 flags to 0.
%
% 2.) Classfication of the components as "good" or "bad".
%     - The classification is done by a user-specified function.
%     - The generic default function is designed for maximum generality.
%     - You can write and use your own classification function. Use
%       identify_artefact_components.m to get input and output formats
%       correct. use S.identfun = @my_custom_function.
%
% 3.) Removal of bad components from the MEG data.
%     - This is a fully automated stage.

for subnum=1:length(spm_files),
        
    % Stage 1 Parameters
    S=[];
    S.to_do = [1 1 1];                              % to_do(1): Do the ICA decomposition. to_do(2): Classify bad components. to_do(3): Remove bad components. 
    S.fname=spm_files{subnum};                      % Input SPM object name (including the path).
    S.logfile = 1;                                  % Set to one to output a log file as well.
    S.overwrite = 0;                                % Set to 1 to overwrite previous ICA runs and re-estimate ICA decomposition.
    S.ica_file = [S.fname '_preproc_ica_results'];  % This the file to which the ICA decompositon is stored.
    S.used_maxfilter = 1;                           % Set to 1 for Maxfiltered Data. Set to zero otherwise.
    
    % Stage 2 parameters
    S.eyetracker_file=eye_tracker_files{subnum};    % The filename of the .txt output from the Eyetracker.
    S.eyetracker_chan=319;                          % The channel in the SPM object that has the Eyeblinks recorded.
    S.do_plots=0;                                   % Set to 1 to produce summary plots of bad components.
    S.manual_approval=1;                            % Set to 1 to activate manual classification of bad components.
    S.do_kurt=4;                                    % Set to -1 to rank all components by kurtosis and visually inspect one by one
                                                    % Set to positive integer to view N (e.g. 5) most and least kurtotic components.
                                                    % Set to 0, otherwise (default).
                                                    
    S.do_mains= 1;                                  % Set to 1 to automatically classify 50Hz mains interference.
    S.do_blinks= 1;                                 % Set to 1 to automatically classify blink artefacts (Eyetracker or EOG required).
    S.do_ecg= 1;                                    % Set to 1 to automatically classfiy ECG artefact (ECG required).
    S.just_ica=0;
    [spm_files_new{subnum}]=osl_africa(S);

end

%% WHAT TO SEE AND DO

% 1) - You will see osl_africa call fastICA and perform the ICA
% decomposition. This will take several minutes so have a careful look at
% the settings we have provided. Think about whether you will use these or
% something different.

% 2) - After the ICA, the result will be saved to disc and then passed to
% identify_artefact_components_manual.m. As we have set S.do_plots = 1 and
% S.manual_approval = 1, a plot of the automatically identified components
% will pop up. Make the fiugre nice and big and look at the different 
% sub-plots.
% 
% You will see the following message:
%
% Please approve artefacts for removal  (e.g. [1 3 5] to select subset, -1 to select all for removal, enter to keep all artefacts):     
%
% If you think that all the components shown are artefacts, enter -1 in the 
% command line and press enter.
%
% Alternatively if you only think the 1st and 3rd components shown are
% artefacts (i.e the top row and 3rd from top row of subplots) then type 
% [1 3] into the command line and press enter.

% 3) - do_kurt - We may not have any signals like ECG to use to identify
% the artefacts. As a result, we have to visually inspect every time
% course. However, we can help ourselves a bit by reordering the
% components. ICA will produce them in a random order but we know that
% blink and heart beat related artefacts tend to have very high kurtosis,
% whilst mains has very low kurtosis. 
% 
% If we set do_kurt to -1, we will visualise each component in order of
% descending kurtosis. Within the first 10 we will probably find the blink
% and heart-beat related components. If we then get bored we can skip the
% remainder.
% 
% Alternatively we can visualise a fixed number of components, by setting
% do_kurt to an integer (e.g. 4). This will show the four lowest kurtosis
% components - which will likely contain a mains component - & then the
% four highest kurtosis components - that will likely have either blink or
% heart beat or both.
%
% For each group, you will see the following message.
%
% Please approve artefact components based on kurtosis for removal  (e.g. [1 3 5] to select subset, -1 to select all for removal, enter to keep all artefacts):  
%
% In the same way as before, identify and flag the component you wish to
% remove.
% 
% Once you have done this, osl_africa will produce a de-noised SPM object.
% Note that if you don;t flag any components as bad then the new file will
% be the same as the old one.
%
% But we have flagged bad components and now we are going to see if this
% has improved our data.
%
%% Let's check the results of AfRICA. Run this cell.
% Load in the the orginal and the denoised data
Dold = spm_eeg_load(spm_files{1});
Dnew = spm_eeg_load(spm_files_new{1});

%% First let's look at the EOG artefact. Run this cell

rho_blinks = zeros(306,1);
for i=1:306;
    rho_blinks(i) = corr(Dold(319,:,1)',Dold(i,:,1)');
end
[max_corr,ind]=max(abs(rho_blinks));

figure; plot(Dold.time,Dold(ind,:,1));
ho; plot(Dnew.time,Dnew(ind,:,1),'r--');
axis tight;
xlabel('Time (s)'); legend('Original Data', 'Data post-AfRICA');
title({'Channel most corrupted by EOG, before and after De-noising' ['Raw channel correlation with Eyetracker ' num2str(max_corr) '.']...
    ['De-noised channel correlation with Eyetracker ' num2str(corr(Dold(319,:,1)',Dnew(ind,:,1)')),'.']});

% Have we reduced the EOG artefact....?

%% Now let's look at the ECG artefact! Run this cell.

rho_ecg = zeros(306,1);
for i=1:306;
    rho_ecg(i) = corr(Dold(find(strcmp(D.chantype,'ECG')),:,1)',Dold(i,:,1)');
end
[max_corr,ind]=max(abs(rho_ecg));

figure; plot(Dold.time,Dold(ind,:,1));
ho; plot(Dnew.time,Dnew(ind,:,1),'r--');
axis tight;
xlabel('Time (s)'); legend('Original Data', 'Data post-AfRICA');
title({'Channel most corrupted by ECG, before and after De-noising' ['Raw channel correlation with ECG ' num2str(max_corr) '.']...
    ['De-noised channel correlation with ECG ' num2str(corr(Dold(find(strcmp(D.chantype,'ECG')),:,1)',Dnew(ind,:,1)')),'.']});

% Have we reduced the ECG artefact....?

%% Finally let's look at the mains interference. Run this cell.

magnetomer_average_orig = mean(Dold(find(strcmp(Dold.chantype,'MEGMAG')),:,1),1);
magnetomer_average_new = mean(Dnew(find(strcmp(Dnew.chantype,'MEGMAG')),:,1),1);

mag_old_fft=fft(magnetomer_average_orig); mag_old_fft = abs(mag_old_fft(1:floor(numel(mag_old_fft)/2)));
mag_new_fft=fft(magnetomer_average_new);  mag_new_fft = abs(mag_new_fft(1:floor(numel(mag_new_fft)/2)));

freq_ax=0:(D.fsample/2)/(numel(mag_old_fft)-1):(D.fsample/2);

figure; plot(freq_ax,mag_old_fft);
ho; plot(freq_ax,mag_new_fft,'r--');
axis tight;
xlabel('Frequency (Hz)'); legend('Original Data', 'Data post-AfRICA');
xlabel('Time (s)'); legend('Original Data', 'Data post-AfRICA');
title({'Spectra of average over all magnetometers, before and after de-noising'});

% Have we reduced the mains artefact....?

%% It is important to note the following things about AfRICA.

% 1.) AfRICA will ignore any bad trials/bad epochs/bad channels. Do not
% change any bad trials/bad epochs/bad channels to good after applying
% AfRICA. You can always set additional trials/epochs/channels to bad
% though.

% 2.) AfRICA uses fastICA which can fail to converge sometimes. If fastICA
% fails to converge, you should not use the de-noised data. Check the log
% file which records the fastICA messages.

% 3.) When you remove bad components from your data, the data will have a
% reduced rank. Any source reconstruction methods will have to account for
% this.
